% function for fitting size of ellipse formed by a tracer-filled diving
% vessel
%
% inputs: 
% im: an image of the vessel
% referenceTiltAngle: the estimate of the tilt angle of the vessel's minor
% axis; if the fit's tilt angle is too far off, then use the reference tilt
% angle
%
% outputs:
% fwhm_minor: estimate of the fwhm of the minor axis of the ellipse
% fwhm: the fwhm through all the angles in te radon transform

function [filtim,bw,fwhm_minor,tiltang,fwhm] = radonFilterDivingVessel_v2_0(im,referenceTiltAngle)

%% radon transform
%generate radon transform and normalize:
radonthresh = .35;
tilttolerance = 8;     %in degrees the allowable deviation from the reference tilt angle
theta = 0:179;
im = im -mean(im(:));
rim = radon(im,theta);
minr = min(rim);
minr = repmat(minr,size(rim,1),1);
maxr = max(rim);
maxr = repmat(maxr,size(rim,1),1);
rnorm = (rim-minr)./(maxr-minr);

%determine fwhm
[~,maxind] = max(rnorm);
[~,Y] = meshgrid(1:size(rnorm,2),1:size(rnorm,1));
temp = repmat(maxind,size(rnorm,1),1);

btm = Y>temp & rnorm<radonthresh;
top = Y<temp & rnorm<radonthresh;

temp = cumsum(btm);
[row,col] = find(temp==1);
[colI,I] = unique(col);
btmrow = nan(length(theta),1);
btmrow(colI) = row(I);

temp = cumsum(top,1,'reverse');
temp = flipud(temp);
[row,col] = find(temp==1);
[colI,I] = unique(col);
toprow = nan(length(theta),1);
toprow(colI) = row(I);
toprow = size(top,1)-toprow;

%% for low signal:noise images, can get some lines that have basically no signal, so can backfill those
toprow(isnan(toprow)) = nanmedian(toprow);
btmrow(isnan(btmrow)) = nanmedian(btmrow);

bwradon = zeros(size(rnorm));
for n = 1:size(rnorm,2)
    cur = zeros(size(rnorm,1),1);
    cur(toprow(n):btmrow(n)) = 1;
    bwradon(:,n) = cur;
end

%% estimate minor axis from radon space:
fwhm = btmrow-toprow;
temp = smooth(fwhm,'sgolay',1);
[~,tiltang] = min(temp);
if nargin > 1 && ~isempty(referenceTiltAngle)
    if abs(tiltang-referenceTiltAngle) > tilttolerance
        fwhm_minor = temp(referenceTiltAngle);
        tiltang = referenceTiltAngle;
    else
        fwhm_minor = prctile(fwhm,10);
    end
else
    fwhm_minor = prctile(fwhm,10);
end
        

% inverse transform:
filtim = iradon(bwradon,theta,'hamming');
bw = binarizeFiltered(filtim);

%trim the images to match the input image
filtim = filtim(1:size(im,1),1:size(im,2));
bw = bw(1:size(im,1),1:size(im,2));

function bw = binarizeFiltered(X)

Xmin = min(X(:));
Xmax = max(X(:));
if isequal(Xmax,Xmin)
    X = 0*X;
else
    X = (X - Xmin) ./ (Xmax - Xmin);
end

% Threshold image - global threshold
bw = imbinarize(X);

return;